% P = BB98d_generate_data(F,dF_dt,P);
% 
% Extended BB98 model that is actually intergrated:
% dT_o/dt = A dT_a/dt - B T_o + C T_a
% Where:
% A = k gamma_a / gamma_o
% B = (  lambda_o + x_o + k x_a) [/ gamma_o]
% C = (k lambda_a + x_o + k x_a) [/ gamma_o]
% Parameters in upper case are consistent with Chan et al (JClim, 2022)
% In this code, however, parameters are in the lower case, 
%               and have the following mapping :: A -> c, B -> a, C -> b

% Note that [--] denotes that the expression is not implimented due to 
% accuracy concerns during optimization.  Rather, gamma_o is fixed to 
% 2e8 J/m^2/K.  Parameters a, b, and c are magnified proportionally.
% 
% Underlying physical model:
% gamma_a dT_a/dt = -lambda_a T_a + x_a (T_o - T_a) + F
% gamma_o dT_o/dt = -lambda_o T_o + x_o (T_a - T_o) + k F
%
% [Reference]
% Barsugli, J. J., & Battisti, D. S. (1998). The basic effects of
% atmosphere?ocean thermal coupling on midlatitude variability. Journal of
% the Atmospheric Sciences, 55(4), 477-493.
% 
% Duo Chan -- Oct. 5, 2022 -- WHOI, Falmouth

function P = BB98d_generate_data(T_a,dT_a_dt,P)

    T_o = zeros(size(T_a,1)+1,size(T_a,2));
    
    if isfield(P,'amp_a')
        do_season = 1; 
    else
        do_season = 0; 
    end

    for ct = 1:size(T_a,1)
        if do_season == 0
            a     = P.a;
            b     = P.b;
            c     = P.c;
        elseif do_season == 1
            f_a   = 1 + P.amp_a * sin(ct/720*2*pi + P.phi_a);
            a     = P.a * f_a;
            f_b   = 1 + P.amp_b * sin(ct/720*2*pi + P.phi_b);
            b     = P.b * f_b;
            f_c   = 1 + P.amp_c * sin(ct/720*2*pi + P.phi_c);
            c     = P.c * f_c;
        end
        dy        =  - a * T_o(ct,:) + b * T_a(ct,:) + c * dT_a_dt(ct,:) * 2e8;
        dy        = dy / 2e8;
        P.dy(ct)  = dy * 2e8;
        T_o(ct+1,:) = T_o(ct,:) + dy * 86400/2;
    end

    N_yr   = size(T_o(2:end,:),1)/2/360;
    N_runs = size(T_o,2);

    P.T_o_m = squeeze(nanmean(reshape(T_o(2:end,:),60,12*N_yr,N_runs),1));
    if size(P.T_o_m,1) == 1, P.T_o_m = P.T_o_m'; end

    P.T_o   = T_o(2:end,:);
    P.t_m   = (1+1/24):1/12:(N_yr+1);    

end

